#include "kernel/feObserver.h"
#include "conductivelayer.h"
#include "EdsCalcTask.h"
#include "LineEdsCalcTask.h"
#include "LineAbDipole.h"
#include "mfAddFuncCalcTask.h"
#include <stdio.h>
#include "Timer.h"

#include "AbDiopleReverseTask.h"
#include "FastMath.h"
#include "VmdReverseTask.h"


void main()
{
	feSimpleObserver obs;
	feObserverManager::GetInstance().Register(&obs);
	


	//check main field's file
	if(!gFileExists(MF_FILENAME_ZDUDT) || !gFileExists(MF_FILENAME_ZU) || !gFileExists(MF_FILENAME_MESH))
	{
		//calc main field
		if(!mfAddFuncCalcTask::CalcMainField())
			return;
	}

	double x = sin(1.0);


	const int n = 40;
	double times[n];
	double step = 1e-6;
	double tau = 1.2;
	times[0] = step;
	for(int i=1; i<n; i++)
	{
		times[i] = times[i-1] + step;
		step *= tau+x;
	}
// 	double times[] = {
// 		
// 	};	
// const int n = sizeof(times)/sizeof(times[0]);
// for(int i=0; i<n; i++)
// {
// 	times[i] /= 1000;	
// }

	ConductiveLayer l[6];


	l[0].H = 20000;
	l[0].Sigma = 1;

	l[1].H = 30;
	l[1].Sigma = 1/100.0;


	l[2].H = 50;
	l[2].Sigma= 1.0/50.0;

	l[3].H = 20;
	l[3].Sigma = 1/10.0;

	l[4].H = 100000;
	l[4].Sigma = 1/50.0;

	l[5].H = 10000;
	l[5].Sigma = 0.01;

{



	LineAbReverseTask t;
	//EdsReverseTask t;

	double resx[n];
	double resy[n];
	double resz[n];


	fePolyline<double> g;
	fePolyline<double> r;
	g.Allocate(10);
	r.Allocate(10);
	double side = 10;

	g.Push(-side/2, -side/2);
	g.Push(side/2, -side/2);
	g.Push(side/2, side/2);
	g.Push(-side/2, side/2);

	double x = 20;
	side = 1;

	r.Push(x-side/2, -side/2);
	r.Push(x+side/2, -side/2);
	r.Push(x+side/2, side/2);
	r.Push(x-side/2, side/2);


	VmdEdsCalcTask *task = new VmdEdsCalcTask();
	//LineEdsCalcTask *task = new LineEdsCalcTask();

	task->Init(l, 1, &g, &r);
	task->Eds(0, times, n, resy);
	FILE *fp = fopen("eds_simple1", "w");
	fprintf(fp, "Master thesis\n\n");
	for(int it=0; it<n; it++)
	{
		fprintf(fp, "%e %e\n", 1000*times[it], 1000*resy[it]);
	}
return;
	LineAbDipole d;	
	d.SetMedium(l, 5);
	
	d.SetLength(1);
	d.Ex(10, 10, 0, resx, times, n);
	d.DBzDt(10, 10, 0, resy, times, n);
	d.Ey(10, 10, 0, resz, times, n);
	ReverseTask::Curve c;
	

	ReverseTask::Curve c1;
	ReverseTask::Curve c2;
	c.Init(0, n);
	c1.Init(0, n);
	c2.Init(0, n);	
	for(int i=0; i<n; i++)
	{
		c[i].X = times[i];
		c[i].Y = resx[i];

 		c1[i].X = times[i];
 		c1[i].Y = resy[i];
 
 		c2[i].X = times[i];
 		c2[i].Y = resz[i];
	}

	ConductiveLayer tl[5];

	tl[0].H = 100;
	tl[0].Sigma = 0.1;

	tl[1].H = 100;
	tl[1].Sigma = 0.1;

	tl[2].H = 100;
	tl[2].Sigma = 0.1;

	tl[3].H = 100;
	tl[3].Sigma = 0.1;

	tl[4].H = 10000;
	tl[4].Sigma = 0.1;

	ReverseTask::Curve cs[] = {c, c1, c2};
	t.BindData(tl, 5, cs, 1, 0, 0);
	
	double eps = 1;
	double eps1 = eps + 10;
	double eps0 = eps1;
	int iter = 0;
	Timer timer;
	timer.Reset();
	while(eps > 1e-2)
	{
		eps = 1;
		eps1 = eps + 10;

		while(fabs(eps1-eps)/max(eps1,eps) > 1e-7)
		{
			t.Iterate();
			eps1 = eps;
			eps = t.GetEps();
			iter++;
			printf("\nIter %d. eps: %e\n", iter, eps);

			printf("Done. Medium:\n");
			int nl=0;
			t.GetSolvation(tl, nl);
			for(int i=0; i<nl; i++)
			{
				printf("%lf %lf\n", tl[i].H, tl[i].Sigma);
			}

			if(eps < 1e-3) break;

		}
		
		double e = 0;
		for(int i=0; i<2; i++)
		{
			e += (tl[i].Sigma - l[i].Sigma)*(tl[i].Sigma - l[i].Sigma);
			if(i!=1)e+=(tl[i].H - l[i].H)*(tl[i].H - l[i].H);
		}
		printf("\nepsilon: %lf\n", sqrt(e));
		
 		if(eps > 1e-2)
		{
			t.SplitMesh();		
			printf("\nSplit! New medium below. Press any key.\n");
			int nl=0;
			t.GetSolvation(tl, nl);
			for(int i=0; i<nl; i++)
			{
				printf("%lf %lf\n", tl[i].H, tl[i].Sigma);
			}
			printf("\n");
		//	getchar();
		}
	}
	printf("\ntime: %lf", timer.Seconds());
	return;
}/*
	l[0].H = 100;
	l[0].Sigma = 1.0/100.0;

	l[1].H = 100;
	l[1].Sigma = 1.0/20.0;

	l[2].H = 1000;
	l[2].Sigma = 1.0/300;

	l[3].H = 500;
	l[3].Sigma = 1.0/10.0;

	l[4].H = 100000;
	l[4].Sigma = 1.0/1000.0;

	l[5].H = 10000;
	l[5].Sigma = 0.01;

	fePolyline<double> g;
	fePolyline<double> r;
	g.Allocate(10);
	r.Allocate(10);
	double side = 500;

	g.Push(-side/2, -side/2);
	g.Push(side/2, -side/2);
	g.Push(side/2, side/2);
	g.Push(-side/2, side/2);

	double x = 0;
	side = 1;

	r.Push(x-side/2, -side/2);
	r.Push(x+side/2, -side/2);
	r.Push(x+side/2, side/2);
	r.Push(x-side/2, side/2);


	//VmdEdsCalcTask *task = new VmdEdsCalcTask();
	LineEdsCalcTask *task = new LineEdsCalcTask();
	double eds[n];
	task->Init(l, 5, &g, &r);
	task->Eds(0, times, n, eds);

	FILE *f = fopen("eds1", "w");

	for(int it=0; it<n; it++)
	{
		fprintf(f, "%e %e\n", 1000*times[it], 1000*eds[it]);
	}
	fclose(f);
return;*/
l[0].H = 10;
l[0].Sigma = 0.1;

l[1].H = 10;
l[1].Sigma = 0.5;

l[2].H = 100;
l[2].Sigma = 0.1;

l[3].H = 100;
l[3].Sigma = 0.5;

l[4].H = 100000;
l[4].Sigma = 0.01;

l[5].H = 10000;
l[5].Sigma = 0.01;
	LineAbDipole dip;
	double resx[n];
	double resy[n];
	double resz[n];
	double xs[] = {10, 20, 30, 40, 50, 60, 70, 80, 90, 500};
	double ys[] = {10, 20, 30, 40, 50, 60, 70, 80, 90, 500};
	dip.SetMedium(l, 5);

	FILE *fp = fopen("out\\out.txt", "w");

	dip.SetLength(200);
	for(int i=9; i<10; i++)
	{
		for(int j=9; j<10; j++)
		{			
			printf("%d %d \n", i, j);
			fprintf(fp, "%d %d\n", int(xs[i]), int(ys[j]));
			fprintf(fp, "%e %e 0\n", xs[i], ys[j]);
			fprintf(fp, "1\n1 1 1 1\n1\nTime(s) Ex(V/M) Ey(V/M) Zero\n");
			Timer t;
			t.Reset();

			dip.Ey(xs[i], ys[j], 0, resy, times, n);
			dip.Ex(xs[i], ys[j], 0, resx, times, n);
			dip.By(xs[i], ys[j], 0, resz, times, n);

			printf("time: %e sec\n", t.Seconds());

			char buf[500];
			sprintf(buf, "out\\ex_%e_%e.txt", xs[i], ys[j]);
			FILE *fex = fopen(buf, "w");

			sprintf(buf, "out\\ey_%e_%e.txt", xs[i], ys[j]);
			FILE *fey = fopen(buf, "w");

			sprintf(buf, "out\\ez_%e_%e.txt", xs[i], ys[j]);
			FILE *fez = fopen(buf, "w");

			for(int it=0; it<n; it++)
			{
				fprintf(fp, "%e %e %e 0\n", times[it], resx[it], resy[it]);
				fprintf(fex, "%e %e\n", times[it], resx[it]);
				fprintf(fey, "%e %e\n", times[it], resy[it]);
				fprintf(fez, "%e %e\n", times[it], resz[it]);
			}

			fclose(fex);
			fclose(fey);
			fclose(fez);
		}
	}


	fclose(fp);
	


	//printf("Done. Medium:\n");
// 	int nl=0;
// 	t.GetSolvation(l, nl);
// 	for(int i=0; i<nl; i++)
// 	{
// 		printf("%e %e\n", l[i].H, l[i].Sigma);
// 	}
	
	feObserverManager::GetInstance().Unregister(&obs);
}